Genomic and Immunological Characterization of Pyroptosis in Lung Adenocarcinoma

Pyroptosis is a programmed cell death that may either promote or hinder cancer growth under different circumstances. Pyroptosis-related genes (PRGs) could be a useful target for cancer therapy, and are uncommon in lung adenocarcinoma (LUAD). The expression profiles, mutation data and clinical information of LUAD patients were included in this study. A pyroptosis-related prognostic risk score (PPRS) model was constructed by performing Cox regression, weighted gene co-expression network analysis (WGCNA), and least absolute shrinkage and selection operator (LASSO) analysis to score LUAD patients. Somatic mutation and copy number variation (CNV), tumor immunity, and sensitivity to immunotherapy/chemotherapy were compared between different PPRS groups. Clinical parameters of LUAD were combined with PPRS to construct a decision tree and nomogram. Red module was highly positively correlated with pyroptosis. Seven genes (FCRLB, COTL1, GNG10, CASP4, DOK1, CCR2, and AQP8) were screened from the red module to construct a PPRS model. Significantly lower overall survival (OS), higher incidence of somatic mutation and CNV, elevated infiltration level of the immune cell together with increased probability of immune escape were observed in LUAD patients with higher PPRS, and were more sensitive to Cisplatin, Docetaxel, and Vinorelbine. We constructed a new PPRS model for patients with LUAD. The model might have clinical significance in the prediction of the prognosis of patients with LUAD and in the efficacy of chemotherapy and immunotherapy.


Introduction
Pyroptosis is a type of cell death programmed caused by the family of proteins known as Gasdermin, which results in cell enlargement, dissolution of plasma membranes, fragmentation of chromosomes, and release of intracellular proinflammatory molecules, thereby triggering inflammation and immune responses [1][2][3]. e relationship between cancer and pyroptosis is a prominent subject in immunology at present. Pyroptosis has a crucial role in enhancing or inhibiting several cancers types, including breast cancer, gastric cancer, esophageal cancer, cervical cancer [4]. In addition, cancer cell pyroptosis can be induced during cancer therapy, including chemotherapy, the treatment by small molecule drugs, and nanodrugs [5]. Recent studies focused on those affecting pyroptotic inflammasomes and promoting pyroptosis molecules, which are expected to be effective targets for the treatment of different cancers [4].
In non-small cell lung cancer (NSCLC) cell lines, simvastatin was found to suppress cancer cell proliferation and migration through inducing pyroptosis [6]. Especially, the gasdermin D (GSDMD) and gasdermin E (GSDME) are two important executioners in the pyroptosis mechanism induced by cancer therapy [7]. A number of pyroptosis core proteins are associated with prognosis of many cancer types such as hepatocellular carcinoma, colorectal cancer, gastric cancer, and lung cancer [3]. Gao et al. discovered that knocking down GSDMD could restrict NSCLC cell growth both in vitro and in vivo, and GSDMD overexpression was significantly associated with poor prognosis in lung adenocarcinoma (LUAD) [8].
To date, several bioinformatics-based studies have identified pyroptosis-related genes (PRGs) in specific cancers. Chen et al. developed a risk model consisting of 6 PRGs, which can successfully be used to evaluate the survival and prognosis of hepatocellular carcinoma and distinguish the risk and predict the immune infiltration and treatment efficiency of HCC [9]. Recent reports provide a novel PRGs signature to predict breast cancer patients' tumor immune microenvironment and prognosis [10]. Zhou's study identified a group of PRGs that can effectively predict ovarian cancer patients' response to chemotherapy and immunotherapy [11]. Luo et al. screened seven possible biomarkers to predict the prognosis of patients with colorectal cancer and provide therapy recommendations for these patients [12]. Lung cancer is a leading cause of cancer-related deaths [13], with lung adenocarcinoma (LUAD) accounting for 40% of the incidence of all lung cancer cases [14]. e role of PRGs depends on the type of cancer, and few PRGs have been found in LUAD. erefore, PRGs play an important role in LUAD.
By performing PRGs bioinformatics analysis, we investigated the LUAD genetic variation on the basis of PRGs in the present research. e pyroptosis-related prognostic risk score (PPRS) model was developed based on the least absolute shrinkage and selection operator (LASSO) and Weighted gene co-expression analysis (WGCNA) regression analysis. e features of differential mutation, biological process, immune infiltration, immunotherapy, and chemotherapy response between PPRS groups were studied. In addition, PPRS was combined with clinicopathological features in the construction of a decision tree and nomogram to optimize the predictive accuracy of the risk of LUAD.

Data Collection and Processing.
e workflow of this study was shown in Figure S1. We obtained the expression profile of gene data, copy number variation (CNV), and somatic mutation data of the treated primary LUAD samples from e Cancer Genome Atlas (TCGA) database (https:// portal.gdc.cancer.gov/). After preprocessing, 500 primary LUAD samples remained in the TCGA cohort. Two other independent LUAD cohorts, GSE31210 and GSE72094, were obtained from the Gene Expression Omnibus (GEO, http:// www.ncbi.nlm.nih.gov/geo/). After preprocessing, 226 samples and 398 samples remained in the GSE31210 and GSE72094 cohorts.

2.3.
e Relation between Pyroptosis Score and LUAD Prognosis. Based on the expression level of 27 PRGs, the pyroptosis score of each sample in LUAD was quantified by ssGSEA and arranged in ascending order. e relation of pyroptosis score, clinical features (including T stage, age, N stage, M stage, sex, survival status, and clinical stage) with overall survival was evaluated by multivariate and univariate Cox regression analysis. Pyroptosis score in different clinical features was analyzed by Wilcoxon test or Kruskal-Wallis test.

WGCNA.
With the R package "WGCNA," a gene coexpression network was developed [17] using gene expression value in the identification of the co-expression gene module. First, the scale-free topology fit index for 1 to 30 powers was computed using the "pickSoft reshold" function. According to blockwiseModules, automatic block module detection was performed. When the independence degree reached 0.8, the appropriate power value was determined, and the module's minimum number of genes was set to 30. Subsequently, highly related modules were merged to form a novel module (parameters: deepSplit � 2, min-ModuleSize � 30, height � 0.25). e highly related modules were merged into a new module (parameters: height � 0.25, deepSplit � 2, minModuleSize � 30). e correlation between eigengene module and PRGs was used to estimate the module-pyroptosis association to identify pyroptosis-related gene modules.

PPRS Was Constructed to Assess the Different Risks of LUAD.
e hub gene of the module was identified by the gene expression of the pyroptosis-related module together with the Pearson correlation analysis on the pyroptosis score. e hub genes not substantially correlated with the survival of patients with LUAD (P > 0.05) were removed by performing univariate Cox regression analysis. e genes showing close correlation with the prognosis of LUAD patients were selected in performing LASSO and multivariate Cox regression analysis. Finally, the screened genes were utilized as variables in constructing the model: pyroptosis-related prognostic risk score (PPRS) � Coefficient(mRNAi) * Expression i. Notably, i refers to the final screened genes.

Functional Enrichment Analysis.
According to the risk score of the median, the samples were separated into two groups, namely the low-risk group and the high-risk group. e R software package "GSVA" was utilized to compute each sample's ssGSEA score in various functions, and the Pearson correlation with PPRS was analyzed for the LUAD samples in TCGA cohort. Furthermore, gene set enrichment 2 Journal of Oncology analysis (GSEA) was performed for two groups based on candidate gene sets in the Hallmark database [15].

Determination of Immune Score and Stromal
Score. e ESTIMATE algorithm [18] was run to estimate immune score, stromal score, and ESTIMATE score in the TME, where ESTIMATE score was the combined score of immune score and stromal score.

Evaluation of Immune Cell Infiltration.
CIBERSORT is a gene expression-based universal deconvolution algorithm that can conduct an estimation of the relative proportion of 22 types of immune cells from tissue gene expression profiles [19]. CIBERSORT (https://cibersort.stanford.edu/) algorithm was used to assess the relative abundance of 22 types of immune cells in 3 LUAD cohorts. Moreover, the differences in the infiltration levels of immune cells between low-risk and high-risk groups were evaluated by performing the Wilcoxon rank-sum test.

Response Analysis regarding Immunotherapy and
Chemotherapy. Based on the data of immune checkpoints obtained from HisgAtlas database [20], the differences between low-risk and high-risk groups in their expression levels got analyzed by performing the Wilcoxon test. e potential clinical effects of immunotherapy in our defined risk group were evaluated by TIDE (http://tide.dfci. harvard.edu/). Additionally, on the basis of the information retrieved from the Genomics of Drug Sensitivity in Cancer (GDSC) database, sample sensitivity to four commonly used clinical chemotherapeutic drugs was predicted. e IC50 determined by R packet "pRRophetic" was used as the comparison index among risk groups.

Establishment of the Decision Tree and Nomogram.
According to the clinical characteristics of LUAD samples in the TCGA cohort such as age, sex, M Stage, N Stage, T stage, Stage and PPRS, decision tree was developed using the R packet "rpart". By Combining these clinical features and PPRS, a nomogram was constructed. e model's performance was then evaluated based on the DCA using the ggDCA package.

Statistical
Analysis. e significant differences in OS between the two PPRS groups were observed by performing the Kaplan-Meier survival analysis and log-rank test. e "time-ROC" software was used to generate ROC. Wilcoxon test was used for the comparison between the two groups. All drawings and statistical analyses were carried out using the R software (R Foundation for Statistical Computing, v.4.0.0). Double-tailed P < 0.05 was considered to be of statistical significance.

PRGs Expression and Genetic Differences in LUAD.
We investigated the presence of somatic mutations in PRGs from e Cancer Genome Atlas (TCGA, on the web: https:// portal.gdc.cancer.gov/)-LUAD samples. Of the 565 LUAD samples from TCGA, 300 PRGs were mutated. e results in the waterfall showed that TP53 had the highest mutation rate (90%), followed by CASP5 (4%), TP63 (4%), CASP1 (3%), and GSDME (3%) (Figure 1(a)). To determine whether these mutations affected the survival and biological function of LUAD, the OS of PRGs wild type and PRGs mutant samples were compared, and we observed no significant change of OS between them (Figure 1(b)). e results of GSEA analysis showed that compared with PRGs wild type samples, E2F targets, MYC targets, G2M checkpoint, mitotic spindle, mTOR signaling, and DNA repair were significantly activated, whereas the p53 pathway was significantly inhibited in PRGs mutant samples (Figure 1(c)). CNV was detected in 27 PRGs, and CNV occurred in all of them, among which the copy number amplification in GSDMD, CHMP6, CHMP4C, CHMP4A, and TP6 was the most obvious, and only copy number loss occurred in CHMP2A and IRF1 (Figure 1(d)).
To study whether the change of copy number had an effect on the expression of 27 PRGs, we analyzed the expression of PRGs in copy number amplification group, copy number deletion group, and copy number no significant change group, and found remarkable differences in the expression levels of these 17 PRGs among three groups (Figure 1(e)). In the 27 PRGs, apart from CASP4, CHMP6, CHMP7, GSDMD, and TP63, the other 22 PRGs were expressed differentially at a significant level between primary tumor and paracancerous samples (Figure 1(f )).

Pyroptosis Is the Main
Factor reatening the Survival of Lung Adenocarcinoma. We ranked the pyroptosis scores of each sample obtained by single sample gene set enrichment analysis (ssGSEA), analyzed the correlation with different clinical features, and found that pyroptosis scores were significantly correlated with N stage and survival state, and that the proportion of dead LUAD samples increased with the increase of pyroptosis scores ( Figure S2A). We performed the multivariate and univariate Cox analyses to explore the influence of each pyroptosis score and clinicopathological characteristic on the prognosis of patients with LUAD. According to the results, it was observed that the prognosis of patients with LUAD was affected by the T stage, pyroptosis score, and N stage in an independent manner ( Figure S2B and S2C). Moreover, the OS duration of the samples with low pyroptosis scores determined based on pyroptosis score grouping was remarkably longer compared to samples experiencing high pyroptosis scores ( Figure S2D). TCGA-LUAD samples were stratified in accordance with several clinical features including Tstage, N stage, survival status, clinical stage, M stage, sex, and age. From the results, it was found that the scores of pyroptosis were different at a significant level among the samples exposed to different survival conditions (either dead or alive), different N stage (N1-N4) and different clinical stage (stage I-IV) stratification, respectively ( Figure S2E).

PPRS Model Construction and Evaluation.
e genes in the red module were filtered. Specifically, the relationship between pyroptosis score and the genes in to the red module was investigated, and the module's hub genes were screened with a criterion of P < 0.01. A total of 73 genes were chosen  for subsequent investigation using the univariate Cox regression analysis (P < 0.05, Figure 3(a) and Table S1). Functional analysis on these screened 73 genes showed that immune-related terms and pathways were significantly enriched, supporting a correlation between pyroptosis and tumor immunity ( Figure S3). Among the 73 chosen genes, 10 of them retained with a minimum λ of 0.0295) by the LASSO-Cox regression model (Figure 3(b)). 7 of the 10 genes were selected to construct the model by stepwise multivariate regression analysis. Among the 7 genes, COTL1, FCRLB, CASP4, and GNG10 were the risk factors, whereas protective factors were CCR2, AQP8 and DOK1 (Figure 3(c)). e PPRS of each LUAD sample in TCGA was calculated and normalized according to our defined formula (PPRS was converted to z-score). LUAD samples were divided into high-PPRS and low-PPRS groups based on the z-score of PPRS � 0. e PPRSs of the samples were arranged in ascending order. e corresponding survival status and the expression of 7 genes showed that the increase of PPRS was accompanied by the increase in the number of dead patients (Figure 3(d)). e high-PPRS group in the TCGA-LUAD cohort had a lower survival rate at the same period, according to the results of the survival analysis (Figure 3(e)).
e areas under the receiver operator characteristic (ROC) curve (AUC) for PPRS were 0.73, 0.7, and 0.66 for the 1-, 3-, and 5-year prognosis prediction, respectively (Figure 3(f )). e same analysis was performed in GSE31210 and GSE72094 cohorts, and the difference of OS between two groups and the efficacy of prediction were validated similarly, thereby confirming the overall accuracy and validity of the PPRS (Figure 3(g), 3(h)).

Genomic Mutation in PPRS Risk
Group. We studied whether PPRS was related to genomic stability, and discovered that some clinical characteristics of patients belonging to the high-PPRS group were more in comparison to those in the low-PPRS group, indicating that the genome of the high-PPRS group was more unstable. is was performed through a comparison of the homologous recombination defects, the number of segments, altered fraction, the score of aneuploidy and tumor mutation burden from both groups (Figure 4(a)). All the five genomic features were positively correlated with PPRS ( Figure 4(b)). e prevalence of CNV together with the somatic mutation in the high-and low-PPRS group was shown in the waterfall map. e prevalence of somatic cell mutation, CNV amplification, and deletion in the high-PPRS group was considerably higher than that in the low-PPRS group (Figure 4(c)).

Enrichment Pathway and Immune Characteristics of PPRS Risk Group.
We determined the function of PPRS and observed that there was a positive correlation between PPRS and cell cycle, replication of DNA, repair of nucleotide excision, homologous recombination, and mismatch repair and other pathways regulating cell proliferation ( Figure 5(a), Table S2). We performed a comparison on the normalized enrichment scores (NESs) present in high-PPRS and low-PPRS groups, and high-PPRS was significantly enriched in MYC targets, E2F targets, G2M checkpoint, mTOR signaling, DNA repair, and other pathways relative to low-PPRS. ese pathways were important pathways that affect cancer cell proliferation ( Figure 5(b)). A total of 22 immune cells types were compared between the high-PPRS group and the low-PPRS group in TCGA-LUAD. We noted that the differences were quite significant in the relative proportion of resting mast cells, monocytes, M0 macrophages, resting dendritic cells, resting memory CD4T cells, resting NK cell, memory B cell as well as activated mast cells between the two groups. e relative proportion of the M0 macrophages and activated mast cells from the high-PPRS group were considerably greater than the low-PPRS group. Moreover, the relative proportion of the other 6 immune cells from the low-PPRS group was found to be higher compared to that from the low-PPRS group ( Figure 5(c)). e Pearson correlation analysis array showed that PPRS was significantly correlated with resting memory CD4Tcells, resting dendritic cells, M0 macrophages, monocytes, and resting mast cells ( Figure 5(d)). Furthermore, immune score, stromal score, and ESTIMATE score were substantially lower in the high-PPRS group in comparison to those in the low-PPRS group ( Figure 5(e)). In addition, the same tumor microenvironment (TME) analysis was also carried out in GSE31210 and GSE72094, and the results can be found in Figure S4.

e Combination of PPRS and Clinicopathological Features Improved the Survival Prediction of LUAD.
e decision tree on the basis of clinical characteristic (M stage, sex, age, N stage Clinical stage, and T stage) together with PPRS showed that only T stage, PPRS, and N stage were retained in the decision tree, and that four different risk subgroups C1-C4 were identified (Figure 7(a)). From C1 to C4, the risk increased gradually, the patients' OS reduced gradually, and a remarkable difference in OS between the groups was observed (Figure 7(b)). Among the four subgroups defined by the decision tree, C1 contained only low-PPRS samples, C2 only included high-PPRS samples, and C3 and C4 samples with high-PPRS accounted for a large proportion (Figure 7(c)). From C1 to C4, the proportion of patients in death status gradually increased (Figure 7(d)). Multivariate and univariate Cox regression analyses of all the PPRS and clinical features indicated that the T stage, PPRS, and N stage were the independent prognostic factors of LUAD (Figure 7(e)).

Discussion
e unclear function of pyroptosis in cancer seems to be contextual and is dependent on genetics, cell type, and the pyroptosis induction duration [21]. e complex effects regarding pyroptosis on the onset and progression of cancer mainly included cancer cell viability inhibition, the influence of cancer cell invasion as well as migration, anti-tumor immunity enhancement, and chemosensitivity enhancement [22]. From the present research, we utilized ssGSEA and univariate and multivariate Cox regression models to identify pyroptosis as the primary risk factor for the overall survival (OS) of patients with LUAD.
Abnormal expression of some important PRGs is often observed in various types of cancer. However, most studies have focused on one or two kinds of PRG, whereas the characterization of the anti-tumor effects is usually a result of the interaction of multiple genes in a highly coordinated manner [11]. In this study, to better quantify the effect of pyroptosis on LUAD, we screened pyroptosisrelated red modules by WGCNA and performed LASSO, univariate Cox regression, and stepwise multivariate Cox regression analysis for the purpose of incorporating seven genes identified in the red module for constructing a PPRS model, which can distinguish the genomic mutation and immune characteristics of patients with different PPRS, and the status of biological pathways. Several genes in the PPRS model have been reported in cancer research.
Coactosin-like protein 1 (COTL1) was reported to be high-expressed in glioma tissues in the study of Shao et al. [23] and is closely correlated with the patient recurrence and prognosis. Functionally, COTL1 enhances the proliferation of cells in vitro and cancer growth in vivo [23]. A study on the peripheral blood mononuclear cells based on peripheral blood RNA-Seq indicated the GNG10 imbalance in the head and neck squamous cell carcinoma, which is related to the survival rate of patients [24]. Wang et al. detected that overexpression of GNG10 promotes the progression of colorectal cancer [25]. CASP4 expresses caspase-4 is a classical regulatory component of pyroptosis [26]. Secretoglobin 3A2-lipopolysaccharide (LPS) can eliminate human colorectal cancer cells by regulating the mechanism of CASP4-related pyroptosis [27]. Shibamoto et al. found that CASP4 expression loss is correlated with the unfavorable prognosis of patients with esophageal squamous cell carcinoma [28]. CCR2 is CC chemokine receptor 2, and CCR2 signal transduction in cancer cells can coordinate the suppression of immune response [29]. AQPs belong to a small membrane transport proteins family, whose abnormal expression plays a role in the onset and progression of several tumors [30], such as in gastric cancer [31], cervical cancer [32], and colorectal cancer [33]. ese studies showed that these genes were tumor markers, and that the coordination between them was likely to have an impact on the development of LUAD.
It is reported that pyroptotic can release tumor antigens and damage-associated molecular patterns, thereby initiating adaptive immunity to enhance the efficiency of immunotherapy [34]. Herein, our analysis results showed that PPRS was not only significantly related to the classical CD274, CTLA4, and PDCD1. And MDSC, T cell exclusion, and TIDE also had significant differences in scores in different PPRS groups. Pyroptotic is also related to    chemotherapy [35]. Recently published studies have shown that Cisplatin induces scorch death through activating the MEG3/NLRP3/caspase-1/GSDMD pathway in triple-negative breast cancer [36]. We observed a positive correlation between the high-PPRS and sensitivity of chemotherapeutic drugs Cisplatin, Docetaxel, and Vinorelbine.

Conclusions
In summary, our study highlighted the importance of pyroptosis in LUAD and observed significantly different expression patterns between normal and LUAD samples. Importantly, we constructed a 7-gene prognostic signature related to pyroptosis, and the signature displayed a favorable performance in predicting LUAD prognosis. Notably, the differences on genomic features, enriched pathways and immune infiltration between PPRS-high and PPRS-low groups demonstrated a potential role of seven prognostic genes in the pyroptosis-related mechanism contributing to LUAD prognosis. e signature offered a comprehensive understanding of the correlation between immunotherapy/ chemotherapy sensitivity of LUAD patients and cell pyroptosis. Our study provides a new insight for understanding pyroptosis-related mechanisms and the hope for developing new therapeutic drugs targeting pyroptosis for LUAD patients.

18
Journal of Oncology